library(numDeriv)

model.selection<-function(data.list,start=c(0,0,0),mode="dynamic") {

base.model.spec<-c(1,1,1,1,1,1,1,1,1,1,1)
base.value<-c(0,0,0,0,0,0,0,0,0,0,0)
# base.value<-c(2,2,2,2,0,0,0,0,0)

par.lower=c(0,0,0,0.0,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf)
par.upper=c(1,1,1,5,Inf,Inf,Inf,Inf,Inf,Inf,Inf)
par.scale=c(1,1,1,1,1,1,1,1,1,1,1)

test.seq<-NULL

#test.seq<-rbind(test.seq,c(0,0,0,0,0,0,0,1,0)) 
# test.seq<-rbind(test.seq,c(0,0,0,0,0,0,0,0,1)) 

test.seq<-rbind(test.seq,c(0,0,0,0,1,0,0,0,0,0,0)) 

test.seq<-rbind(test.seq,c(0,0,0,0,0,1,0,0,0,0,0)) 
test.seq<-rbind(test.seq,c(0,0,0,0,0,0,1,0,0,0,0)) 
test.seq<-rbind(test.seq,c(0,0,0,0,0,0,0,1,0,0,0)) 
test.seq<-rbind(test.seq,c(0,0,0,0,0,0,0,0,1,0,0)) 

test.seq<-rbind(test.seq,c(0,0,0,0,0,0,0,0,0,1,0)) 
test.seq<-rbind(test.seq,c(0,0,0,0,0,0,0,0,0,0,1)) 

test.seq<-rbind(test.seq,c(1,0,0,0,0,0,0,0,0,0,0)) 
test.seq<-rbind(test.seq,c(0,1,0,0,0,0,0,0,0,0,0))
test.seq<-rbind(test.seq,c(0,0,1,0,0,0,0,0,0,0,0)) 
test.seq<-rbind(test.seq,c(0,0,0,1,0,0,0,0,0,0,0)) 

test.ans.list<-list()
current.model.spec<-base.model.spec

st<-start[current.model.spec==1]
ll<-par.lower[current.model.spec==1]
ul<-par.upper[current.model.spec==1]
sc<-par.scale[current.model.spec==1]

current.ans<-fit.model(st,lower=ll,upper=ul,scale=sc,par.restrict=base.model.spec,par.restrict.value=base.value,
                       data.list=data.list)

msg<-"no issue"
if((current.ans$message!="relative convergence (4)")&&(current.ans$message!="X-convergence (3)")&&
   (current.ans$message!="both X-convergence and relative convergence (5)")) msg<-"no convergence 0"

pv.list<-NULL
for(i in 1:nrow(test.seq)) {
  test.model.spec<-current.model.spec-test.seq[i,]*(current.model.spec==1)
	cat(i," ",test.seq[i,]," current model=",current.model.spec,"test mode=",test.model.spec,"\n")
  
# browser()
	start.use<-start
#  start.use<-convert.par.stl(current.ans$par,model.spec=current.model.spec)
#  start.use<-as.numeric(ifelse(start.use=="N/A",start,start.use))
  st<-start.use[test.model.spec==1]

# browser()

  ll<-par.lower[test.model.spec==1]
  ul<-par.upper[test.model.spec==1]
  sc<-par.scale[test.model.spec==1]

#  cat("st=",st,"\n")
  test.ans<-fit.model(st,lower=ll,upper=ul,scale=sc,par.restrict=test.model.spec,par.restrict.value=base.value,
                      data.list=data.list)

	if((test.ans$message!="relative convergence (4)")&&(test.ans$message!="X-convergence (3)")&&
	     (test.ans$message!="both X-convergence and relative convergence (5)")) msg<-paste("no convergence ",i,sep='')

	llk.current<-current.ans$objective
	llk.test<-test.ans$objective
	lrt<-2*(llk.test-llk.current)
  df<-min(1,abs(sum(current.model.spec)-sum(test.model.spec)))
	p.value<-1-pchisq(lrt,df)
  pv.list<-c(pv.list,p.value)
#  if(i==nrow(test.seq))  browser()
  
	cat(test.seq[i,]," current ",llk.current," test ",llk.test," p-value ",p.value,"\n")
	cat(current.ans$message," ",test.ans$message,"\n\n")

test.ans.list[[i]]<-list(current.model.spec,current.ans,test.model.spec,test.ans)
if(p.value>0.05)  if(mode=="dynamic") {
		current.ans<-test.ans
		current.model.spec<-test.model.spec
	}
}

#hm<-hessian(llk,current.ans$par,par.restrict=current.model.spec,par.restrict.value=base.value,data=data)

#if(det(hm)<=1e-5) se<-rep("singular",length(current.ans$par))
#else se<-sqrt(diag(solve(hm)))

se<-rep("N/A",length(current.ans$par))

list(current.model.spec,current.ans,msg,test.ans.list,pv.list,test.seq,base.value,se)
}

convert.par.stl<-function(par,model.spec=NULL) {

	long.par<-rep(-104.3,length(model.spec))
#	long.par[1]<-par[1]	

  counter<-1
  for(i in 1:length(model.spec)) {
    if(model.spec[i]==1) {
      long.par[i]<-par[counter]
      counter<-counter+1
    }
  }
  long.par[long.par==-104.3]<-"N/A"
	long.par
}

summary.stat<-function(table) {
  total.cp<-sum(table[,"model selection msg"]=="no issue")
  a<-subset(table,table[,"model selection msg"]=="no issue")
  model.feature.count<-NULL
  for(i in 5:11) model.feature.count<-c(model.feature.count,sum(a[,i]==1))
  
  label<-c("total",dimnames(table)[[2]][5:11])
  result<-cbind(label,c(total.cp,model.feature.count),c(1,model.feature.count/total.cp))
  
  model.complexity<-apply(matrix(as.numeric(a[,5:11]),ncol=7),1,sum)
  
  result<-list(result,cbind(as.numeric(a[,3]),model.complexity))
  aa<-result[[2]]
  
  aa.low<-subset(aa[,2],aa[,1]<=median(aa[,1]))
  aa.high<-subset(aa[,2],aa[,1]>median(aa[,1]))  
  
  result[[3]]<-wilcox.test(aa.low,aa.high)
  result
}

coef.summary<-function(model.summary) {
  cn<-dimnames(model.summary)[[2]]
  coef.list<-cn[grep("coef",cn)]
  cs<-NULL
  cs2<-NULL
  for(coef in coef.list) {
    a<-as.numeric(model.summary[,coef])
    na.count<-sum(!is.na(a))
    pos.count<-sum(na.omit(a)>0)
    neg.count<-sum(na.omit(a)<0)
    total<-length(a)
    cs<-rbind(cs,c(coef,na.count,pos.count,neg.count,total)) 
    a<-na.omit(a)
    a<-subset(a,(a>quantile(a,probs=0.1))&(a<quantile(a,probs=0.9)))
    cs2<-rbind(cs2,c(coef,summary(a),sd(a)))
  }
  dimnames(cs)[[2]]<-c("coef name","significant","postive","negative","total")
  list(cs,cs2)
}

ms.summary<-function(ms,coef.name=NULL) {
  long.par<-convert.par.stl(ms[[2]]$par,model.spec=ms[[1]])
  
  pv.list<-ms[[5]]
  test.seq<-ms[[6]]
  base.value<-ms[[7]]
  
  result<-cbind.data.frame(coef.name,long.par)
  
  p.value<-rep("NA",nrow(result))
  null.hypothesis<-rep("NA",nrow(result))
  msg1<-rep("NA",nrow(result))
  msg2<-rep("NA",nrow(result))
  
  #  browser()
  
  for(i in 1:nrow(test.seq)) {
    if(sum(test.seq[i,])==1) {
      p.value[test.seq[i,]==1]<-pv.list[i]
      null.hypothesis[test.seq[i,]==1]<-base.value[test.seq[i,]==1]
      msg1[test.seq[i,]==1]<-ms[[4]][[i]][[2]]$message
      msg2[test.seq[i,]==1]<-ms[[4]][[i]][[4]]$message
      
      #      browser()
    } else {
      p.value[test.seq[i,]==1]<-pv.list[i]
      null.hypothesis[test.seq[i,]==1]<-paste(base.value[test.seq[i,]==1]," (test ",i,")",sep="")
      
      msg1[test.seq[i,]==1]<-ms[[4]][[i]][[2]]$message
      msg2[test.seq[i,]==1]<-ms[[4]][[i]][[4]]$message
      
      #            browser()
    }
    #  browser()
    
  }
  
  se<-convert.par.stl(ms[[8]],model.spec=ms[[1]])
  result<-cbind.data.frame(result,se,null.hypothesis,p.value,msg1,msg2,rep(ms[[2]]$objective,nrow(result)))
  dimnames(result)[[2]]<-c("coef","estimate","std err","null (if applicable)","p-value","current msg","test msg","llk")
  
  result
  
}

